The Solution of the Modified Helmholtz Equation in a Wedge 
and an Application to Diffusion-Limited Coalescence 
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Abstract: The general solution of the modified Helmholtz equation, q X x(x,y) + q yy (x, y) 
A(3 2 q(x, y) = 0, in the wedge < x < y < oo, is presented. This solution is used to find 
the explicit steady-state of diffusion-limited coalescence, A + A # A, on the half-line, with 
a trap-source at the origin. 
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A new method for solving boundary value problems for linear and for integrable 
nonlinear PDEs was introduced in [1]. This method has been applied so far to evolution 
equations [2-4] and to Laplace's equation [5]. In this letter we apply this method to the 
equation 

E xx + E yy + 1 {-E x + E y ) = 0, (1) 

in the wedge < x < y, where E(x,y) is a scalar function and 7 is a positive constant. 
Using the substitution E(x,y) = e~^^ v ~ xS) q{x,y) + e~ l( - M ~ x \ eq. (1) becomes the modified 
Helmholtz equation 

q xx + q y y-4f3 2 q = 0, 0<x<y, f3 = -J= . (2) 

Eq. (1) with 7 = v/D represents the steady state of the diffusion-limited reaction A + 
A ^ A on the line, where the A-particles diffuse with diffusion constant D, they merge 
immediately upon encounter, and split into two particles (the back reaction) at rate v [6- 
8]. E(x,y) represents the probability that the interval (x,y) is empty. The concentration 
profile of the particles is related to E through c(x) = —E y (x,x) = 7 — q y (x,x). Suppose 
that we limit ourselves to the positive x-axis, then the domain of Eq. (1) is < x < y. The 
forward reaction is described by the boundary condition (BC) E(x,x) = 1. If the origin 
absorbs all impinging particles, but also generates particles at rate v' , then E x (0,y) = 
7'£'(0,y), where 7' = v'/D. The associated BCs for q(x,y) are 

q(x, x) = 0, qx (Q,y) + (l- 7 ')g(0, y) = (7' - l)e~^ y . (3) 

Particular solutions are known for the case of a perfect trap at the origin (7' = 0), and 
for the case of 7' = 7. In the latter case the steady state is identical to that of an infinite 
system, without a trap. For 7' > 7, it is not clear a priori whether a steady state even 
exists. The present method shows that it does, and it provides an explicit solution for all 

7 > 0. 



The implementation of the method of [1] to eq. (2) involves the following steps: 
(a) Construct a Lax pair for eq. (2). This means finding two linear eigenvalue equations 
whose compatibility condition is eq. (2). (b) Perform the simultaneous spectral analysis 
of this Lax pair. This means constructing a function fi,(x,y, k) which satisfies both equa- 
tions defining the Lax pair, and which is bounded for all complex values of the spectral 
parameter k. This can be used to: (i) find an integral representation of q(x,y) in terms 
of certain functions q a (k) and qb(k), which we refer to as spectral data; (u) find a global 
relation satisfied by the boundary values of q. (c) Analyze the global relation. This can 
be used to express q a (k) and qb(k) in terms of the given data on the boundary. 

Steps (a) and (b) are algorithmic. The Lax pair, the integral representation for q(x, y), 
and the global relation, are given by eqs. (6), (14), (13). We emphasize that these steps are 
implemented without specifying any boundary conditions for eq. (2). The implementation 
of step (c) depends on the particular type of boundary conditions. 

We begin by finding a Lax pair. Let z = x + iy. Then d z = (1/2)(<9 X — id y ), dj = 
(1/2) (d x + id y ). Using these formulas, eq. (2) can be rewritten in the form q z j — (3 2 q = 0. 
This equation is the compatibility condition of the following equations: 

fi z + ikfi = q z , fi zJ - (3 2 fi = . (4) 

Indeed, applying the operator d 2 - — (3 2 to the first of eqs. (4) and using the fact that it 
commutes with d z + ik, we find (d 2 - — (3 2 )q z = (d z + ik)(d 2 - — (3 2 )fx = 0. Differentiating 
the first of eqs. (4) and using the remaining eq. and q z j — (3 2 q = 0, it follows that 

H z -t—H = -i—q . (5) 

Finally, reverting to the x and y variables the first of eqs. (4) and eq. (5) yield 

/j x + t{k- —)/j, = q z -i—q , fi y - (k + —)n = iq z - —q . (6) 

In summary, if q satisfies eq. (2), eqs. (6) are compatible. 
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In what follows we will construct a function fj,(x, y, k) which solves both eqs. (6) and 
which is sectionally analytic in the entire complex plane, including infinity. The function \x 
has different representations in different sectors of the complex /c-plane, namely /io(x,y, k), 
k-R > 0; fj,\(x,y,k), /cr < and k\ > — /cr; fi2(x,y,k), /cr < and k\ < — /cr, see figure 1. 
The functions /j,q, /j,i, \i 2 are defined as follows: 

Ho(x,y,k) = - Ei(y-y',k)[iq z (x,y') - —q(x,y')]dy' , (7) 



f y P 2 

Hj(x,y,k) = E 1 (y-x,k)F j (x,k) + / E x {y -y ,k)[iq z (x,y ) - — q(x,y)]dy, 



J = 1,2, 

(8) 



(9) 



where 

f°° (3 2 
F x (x, k) — — (1 + i) I E 2 (x — x' , k)[q z (x , x) — q(x ,x )] dx , 

Jx h 

El (y, k) = e^^ , E 2 {x, k) = e [(i-0*+(i+0£]* , 

and F 2 is defined by a similar expression with j Q instead of — J . 

In order to derive these formulas we begin with the second of eqs. (6). A particular 

solution of this equation is given by eq. (7). Assuming that q(x,y), q x (x,y) and q y (x,y) 

tend to zero asi-> oo, it follows that /io also satisfies the first of eqs. (6). We note that 

Re{/c + (3 2 /k} = /cr(1 + /3 2 /|/c| 2 ), thus, since y — y' < 0, the exponential appearing in (7) 

is bounded at infinity for /cr > 0. Also Re{,3 2 /&;} = /3 2 /cr/|/c| 2 , thus, since y — y' < 0, this 

exponential is bounded at zero for /cr > 0. The general solution of the second of eqs. (6) is 

given by the r.h.s. of eq. (8), where Fj(x, k) is replaced by F(x, k) = fj,(x,x, k). We choose 

the function F(x, k) in such a way that \i also satisfies the first of eqs. (6): 

dF_ = d^x,k) + d»(x,x,k) = [(1 _ z)fc+(1+z) ^ ]jP+(1+z)fe(x ^ ) _^ (x?x)] . (10) 

The integral term of eq. (8) is bounded for /cr < 0, thus fi is bounded for /cr < provided 
that there exist solutions of eq. (10) which are bounded for /cr < 0. A particular solution 
of eq. (10) is given by eq. (9). We note that the real part of the exponent of E 2 involves 

(k R + ki)(l + /3 2 /\k\ 2 )x, and Re{(l + i) f3 2 / k} = (k R + h) f3 2 / \k\ 2 , thus, since x-x' < 0, the 
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exponential appearing in Fi(x, k) is bounded both at infinity and at zero, provided that 
kn + kj > 0. Hence fj,i(x, y, k) is bounded for /cr < and k\ > — /cr. Similar considerations 
apply to the function Hi- 

The function fx(x,y,k) is a sectionally analytic function of k which satisfies \i = 
—i\ + O(p-), k — »■ oo. Thus, this function can be reconstructed from knowledge of its 
"jumps" across the contour depicted in figure 1. If we use the convention that the © 
region is to the left of the orientation of the contour, these jumps are Hi — Hu A*2 _ A*o> 
and fio — Hi. An important result of the method of [1] is that these jumps, denoted by 
J(x, y, k), have an explicit x and y dependence. Indeed, eqs. (6) imply J x +i(k—(3 2 /k)J = 0, 
J y -(k + [3 2 /k)J = 0. Thus, J(x,y,k) = exp[(jfe + [3 2 /k)y - i(k - [3 2 /k)x]p(k), where p(k) 
is a function of k only and can be computed by evaluating J on the boundary of the given 
domain. This yields 

li 2 -Hi = E 3 q a (k) , -h = k R <0; E 3 (x,y, k) = exp[(k + (3 2 /k)y - i(k - (3 2 /k)x] ; 

A*2 -Mo = E 3 q b (k) , 00 > ki > ; /j ~ Hi = E 3 {q a (k) - q b (k)] , 00 < fci < , 

(11) 

q a (k) = (l + i) / £ , 2 (-a;,A;)[q' 2 (a;,a;) - — q(x,x)]dx , 

° (12) 

/ 8 

Ei(-y,k)[iq z (0,y)- —q(0,y)]dy, 

and Ei(x, k), E 2 (y, k) are defined in (9). The function Hq{ x i x i k) in addition to satisfying 
eq. (7) evaluated at y — x, also satisfies eq. (10), i.e., it is given by 

f°° 1 ft 1 

Ho(x,x, k) — — (1+i) / E2(x—x',k)[q z (x',x') — — q(x',x')] dx' , /cr > and k\ > — /cr . 

Jx k 

Evaluating this equation at x = and using eq. (7) to compute /u(0, 0, k) we find 

Qa(k) = q b (k) , k R > and fci > -/c R . (13) 

This condition implies that \i\ provides the analytic continuation of \iQi ^- e -i there is no 
jump between Ho an d Hi- 
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Eqs. (11), together with the behavior of// for large fe, define a scalar Riemann-Hilbert 
(RH) problem [9,10] whose unique solution is 



loo 



,3iTT/4 



1 /° E 3 (x,y,e)W) je 1 f e °° E 3 (x,y,£)q a (£) J§ 

p{x,y,k) = — — r d£+— — d£ 

2m l_ ir ~ £ — k 2m L £ — k 



where the exponential E 3 (x, y, fe) is defined in eqs. (11). Using eq. (5) to express q(x,y) in 
terms of \x and \x-^ it follows that 

q(x, y) = —f° e <*+£>«-*(*-£)* &M dk + J- f 3W4 °° e (*+£)»-*<*-£>* ^ dk . 

2m J_ ioo k 2m J k 

(14) 
In summary, suppose that q(x, y) satisfies eq. (2). Assume that q and its derivatives decay 
to zero at infinity. Furthermore, assume that appropriate BCs are prescribed at x = and 
at y = x so that there exists a global solution. This global solution is given by eq. (14) 
where q a (k) and <#>(fe) are defined by eqs. (12). Furthermore, the boundary values of q 
satisfy the global relation (13). 

We will now use the above formulation to solve a concrete boundary value problem. 
Let q(x,y) satisfy the BCs 

q(x,x) = f(x), q x (0, y) + aq(0, y) = g(y) , (15) 

where a is a given constant and f(x), g(y) are given functions decaying to zero at large 
x and large y respectively. The expression for q a (k) involves q x (x,x), q y (x,x) and q(x,x). 
Using the first BC in (15), it follows that q(x, x) = f(x), q x (x, x) = f'(x) — q y (x, x). Thus 
q a (k) involves the unknown function q y (x, x) and a known function expressed in terms of 
f'(x) and f(x). Using integration by parts to simplify the latter expression we find 

1 + i (3 2 f°° f°° 

Qa(k) — /(0) + (fe - — ) / E 1 (-x,k)f(x)dx-i E 1 (-x,k)q y (x,x) dx . (16) 

^ k Jo Jo 

Similarly, the expression for qj, involves q x (0,y), Qy(0,y) and q(0,y). Using the second 
BC in (15), it follows that q x (0,y) = g(y) — aq(0,y). Thus <#,(fe) involves the unknown 
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functions q(0,y) and q y (0,y), and a known function expressed in terms of g{y). Using 
integration by parts to simplify the former expression we find 

1 i f°° 1 3 2 Z" 00 

Qb(k) = --f(0) + - / E 2 (-y,k)g(y)dy + -(k-^--ia) / E 2 (-y, %(0, y) dy . (17) 



Hence, eq. (13) becomes 

1 B 2 

-(k - y - ia)Q b (k) + iQ a (k) = F(k) , k K > and h > -k R , (18) 

where Q a (k) is the last integral term in eq. (16), Qb(k) is the last integral term in eq. (17) 
and F(k) is a known function: 



F(k) = -g/CO) + (k - tL) J El (-x, k)f(x) dx--J Eh 



(-y,k)g(y)dy. (19) 



Eq. (18) is one equation for two unknown functions. However, using the analytic properties 
of Q a and Qb it is possible to obtain both these functions through the solution of a scalar RH 
problem. Indeed, we first note that Qb(k) remains invariant under complex conjugation; 
thus, taking the complex conjugate of eq. (18) we find an equation involving Q a {k) and 



Q a (k), valid for \k\\ < /cr, /cr > 0. The exponentials of Q a (k) and of Q a (k) involve 
(1 — i)k + (1 + i)jr and (1 + i)k + (1 — i)jr- Under the transformation k — ► (1 + i)k 

rt2 rt2 

these terms become 2k + %■ and i(2k — j-). Thus, using k — ► (1 + i)k in the equation 



involving Q a (k) and Q a (k) and then following complex conjugation we obtain an equation 

involving J °° q y (x, x) exp[z(2/c — ^-)x] da; and the complex conjugate of this expression, 

valid for /cr > 0. Finally, letting k — > — fe in this equation we find an equation valid for all 

kn: 

A(k)A(-k)- 



where 



$+(fc)- ^ vv ^ v 7 $+(fc)=g(fc), fcr = 0, (20) 

A(-fc)A(fc) 



$+(A;)= f^ e^ 2k -^ x q y {x,x)dx, A{k) = k 2 - (1 + ' )a fc + ^ , (21) 

Jo 2 2 



and B(k) = G{k) for A; R > 0, B{k) = G(-k) for k R < 0, with 



A(-k), A(k) — -=- A(k) 



G(k)=i^^[-^-F((l-i)k) + =^F((l-i)k)-F((l + i)k)-F((l + i)k)}. (22) 
A(k) A (~ k ) A(-k) 

Since Re{i{2k - [3 2 /k)} = ~h(2 + P 2 /\k\ 2 ), Re{-i[5 2 /k} = -ki[3 2 /\k\ 2 , $+(fc) is 



analytic for k\ > 0, thus $+(/c) is analytic for k\ < 0. These facts, together with the 
estimate $ + (/c) = O(^), k — > oo, imply that eq. (20) defines an inhomogeneous scalar RH 
problem. The solution of this problem depends on the index of the jump function. Let 
C + (C~) denote the upper (lower) half of the complex /c-plane. In our case the index is 



given by (number of zeros of A(k) in C + ) + (number of zeros of A(—k) in C + ) — (number 
of zeros of A{— k) in C + ) — (number of zeros of A(k) in C + ). If N^ denotes the number 
of zeros in C 1 * 1 , the index equals 2N + — 2N~ . Thus there exist 3 cases: (a) A(k) has 1 
zero in C + , then the index is zero and eq. (20) has a unique solution, (b) A(k) has 2 
zeros in C + , then the index is 4 and the solution of (20) depends on 3 arbitrary constants. 
(c) A(k) has 2 zeros in C~, then the index is —4 and the solution of (20) exists only if 
the function B(k) satisfies 4 orthogonality conditions. The roots of A(k) = are given by 



4A = (1 + i)a(l ± yj\ - (^) 2 ). Thus if ||| < ^, then iV+ = 2 if a > and N~ = 2 if 
a < 0. If |£| > ^ then A^+ = 1. 

In summary, if | — | > ^Q- then eq. (2) with the BCs (15) has a unique solution; if 



I — I < ^r and a > 0, there exist 3 arbitrary constants; if I — I < ^ and a < there exists 

\ a ' 2 J ' i a i 2 

a solution only if f{x) and g(y) are such that they satisfy 4 orthogonality conditions. 

We now consider the case of iV + = iV - = 1. Let A + and A - denote the roots of A(k) 

in C + and C~ respectively. Then the homogeneous version of eq. (20) becomes 

x+w _(fc-A+)(*-A-)(* + g)(W 
V ; (ife + A+)(A; + A-)(A;-A+)(ife-A-) 

The canonical solution of this problem, i.e., the solution satisfying X(k) = 1 + ^(r)) 
fe — > oo, is given by 



Thus, 



gm - X ^ r m di (25) 



When N = 2, the solution is still given by eq. (25), but now X + (k) = A(k)A(—k). 
In this case, however, the solution exists only if the orthogonality relations 



_^_^ = , n = 0,l,2,3, (26) 

are satisfied. 

In the practical physical application of diffusion-limited coalescence, the BCs (3) imply 
a = 7'/2 — 7, f(x) = 0, g(y) = (7' — 7)e~ ^ y . The three cases for the zeros of A(k) 
discussed above correspond to: (i) 7' < 0, for A^ + = 2, iV~ = 0, (u) < 7' < 7, for 
iV" 1 " = A^~ = 1, and (Hi) 7' > 7, for A^ + = 0, iV - = 2. Thus, case (i) is unphysical. The 
existence of a unique solution for case (ii) is consistent with the fact that at the endpoints, 
7' = 0,7, solutions were already known [7,8]. For case (in), it can be easily verified that 
the orthogonality conditions (26) are satisfied, and hence there exists a unique steady- 
state solution even when 7' > 7. As a final remark, we note that q y (x,x), and hence 
the concentration profile c(x) = 7 — q y (x,x), are already derivable from $(/c). Indeed, 
the substitution i(2k — (3 2 /k) = —s in 3>(h) yields the Laplace transform of q y (x, x), 
j e~ sx q y (x,x) dx, which is then inverted to find the concentration. 
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FIGURE CAPTIONS 

Figure 1: The regions of analyticity of (j,q, fix, \ii in the complex /c-plane. The vertical 
line is /cr = 0, the other ray is k\ = — /cr and /cr < 0. 
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